/***************************************************************************************************************************/
// But de programme est de faire le simulateur d'epidemie SEIR en utilisant la méthode directe de Gillespie en 1977
// Etape1: choisir le temps ou la reaction se produit avec la probabilite p1
// Etape2: choisir une ville ou la reaction se produit avec la probabilite p2
// Etape3: choisir une reaction des 8 reactions produite avec la ville choisie dans l'etape 2 avec la probabilite p3
// Par : TRAN Thi Cam Giang, P15 -IFI
/***************************************************************************************************************************/


/*
 *  BatchVacc.h
 *  BatchVacc
 *
 *  Created by TRAN Thi Cam Giang on 24/11/11.
 *  Copyright 2012 __MyCompanyName__. All rights reserved.
 *
 */

//Initialiser les bibliotheques utilisees
#include "batchVac.h"
#include <ctime>
#include <cstdlib>
#include <stdio.h>
#include <cstring>
#include <cmath>
//#include <algorithm>
#include <vector>
//#include <list>
//#include <map>
//#include <stdlib.h>
//#include <math.h>
//#include "stdio.h"
//#include <stdio.h>
#include <fstream>
#include <iostream>
#include <fstream>
#include <sstream>
//#include<sys/stat.h>
//#include<sys/types.h>
//#include <time.h>
//#include <stdlib.h>

using namespace std;

/******************************************************************************************************/
//variables globales
double tauxVaccToutesVilles = 0.0; // entre 0 et 100
int nombreMaxVaccins = 0; // nombre max de vaccins utilises
int nombreTotalVaccine = 0; //somme de vaccins utilise
int nombreTotalInfectes = 0; //somme de infectes
// La valeur du nombre de Pi
 double Pi = 3.141592654 ;
/******************************************************************************************************/
//Initialiser les valeurs initiales pour les parametres de SIM
vector<double> initialervalParSIM(int argc, char** argv){
    vector<double> res;
    //
    double tmax = 5*365;
    res.push_back(tmax);// tmax = 2000(jours)
    double nbVille = 1;
    res.push_back(nbVille); // nbVill = 1

    /*
    double rinitN = 10000;
    res.push_back(rinitN);
    double rinitI = 100;
    res.push_back(rinitI);

    double beta0 = (double)1000/365;
    res.push_back(beta0);
    */

    double sigma = (double)1/7;
    res.push_back(sigma);

    double gamma = (double)1/7;
    res.push_back(gamma);

    double rmu = (double)1/(10*365);
    res.push_back(rmu);

    double epsilon = 0.1;
    res.push_back(epsilon);

    double topology = 0;
    res.push_back(topology);

    double rho = 0.01;
    res.push_back(rho);

    double unitTemps = 1.0;
    res.push_back(unitTemps);

    double graine = 10;
    res.push_back(graine);
    //  iS0 = 13, iI0 = 14, iE0 = 15, iR0 = 16, iN0 = 17,
    //ibeta0 = beta0, ibeta1 = 19, iphi0 = 20;
    double S0 = 97000;   double E0 = 200; double I0 = 100; double R0 = 0.0; double N0 = 10000;
    res.push_back(S0); res.push_back(E0); res.push_back(I0); res.push_back(R0); res.push_back(N0);

    double beta0 = (double)1000/365;
    res.push_back(beta0);

    double beta1 = 0.25;
    res.push_back(beta1);

    double phiMIN = 0.0;
    res.push_back(phiMIN);
    double phiMAX = (double)Pi/2;
    res.push_back(phiMAX);

   //Parameters pour conserver les entrees
    char str[255] = "";
   for (int i=1; i<argc; i++)
       {
           // Recuperer le premier point de la fonction de transformation
             strcpy(str,argv[i]);
             if ((strcmp(str,"-sigma") == 0))
             {
                 i++;
                 res[isigma] =(double) atof(argv[i]);
             }
             else if((strcmp(str,"-gamma")) == 0)
             {
                 i++;
                 res[igamma] = atof(argv[i]);
             }
             /*
             else if((strcmp(str,"-beta"))==0)
             {
                 i++;
                 res[ibeta] = atof(argv[i]);
             }
             */
             else if((strcmp(str,"-rmu"))==0)
             {
                 i++;
                 res[irmu] = atof(argv[i]);
             }
             /*
             else if((strcmp(str,"-initI"))==0)
             {
                 i++;
                 res[irinitI] = atof(argv[i]);
             }
             else if((strcmp(str,"-initN"))==0)
             {
                 i++;
                 res[irinitN] = atof(argv[i]);
             }
             */
             else if((strcmp(str, "-nbVilles"))==0)
             {
                 i++;
                 res[inbVille]= atof(argv[i]);
             }
             else if((strcmp(str,"-topology"))==0)
             {
                 i++;
                 res[itopology] = atof(argv[i]);
             }
             else if((strcmp(str,"-rho"))==0){
                 i++;
                 res[irho] = atof(argv[i]);
             }
             else if((strcmp(str,"-epsilon"))==0){
                 i++;
                 res[iepsilon] = atof(argv[i]);
             }
             else if((strcmp(str,"-tmax"))==0){
                 i++;
                 res[itmax] = atof(argv[i]);
             }
             else if((strcmp(str,"-unitTemps"))==0){
                 i++;
                 res[iunitTemps] = atof(argv[i]);

             }
             else if((strcmp(str,"-graine"))==0){
                 i++;
                 res[igraine] = atof(argv[i]);
             }
             //  iS0 = 13, iI0 = 14, iE0 = 15, iR0 = 16, iN0 = 17,
             //ib0 = 18, ib1 = 19, iphi0 = 20;
             else if((strcmp(str,"-S0"))==0){
                 i++;
                 res[iS0] = atof(argv[i]);
             }
             else if((strcmp(str,"-E0"))==0){
                 i++;
                 res[iE0] = atof(argv[i]);
             }
             else if((strcmp(str,"-I0"))==0){
                 i++;
                 res[iI0] = atof(argv[i]);
             }
             else if((strcmp(str,"-R0"))==0){
                 i++;
                 res[iR0] = atof(argv[i]);
             }
             else if((strcmp(str,"-N0"))==0){
                 i++;
                 res[iN0] = atof(argv[i]);
             }
             else if((strcmp(str,"-beta0"))==0){
                 i++;
                 res[ibeta0] = atof(argv[i]);
             }
             else if((strcmp(str,"-beta1"))==0){
                 i++;
                 res[ibeta1] = atof(argv[i]);
             }
             else if((strcmp(str,"-phiMIN"))==0){
                 i++;
                 res[iphiMIN] = atof(argv[i]);
             }
             else if((strcmp(str,"-phiMAX"))==0){
                 i++;
                 res[iphiMAX] = atof(argv[i]);
             }

        }

return res;

}

void coutEcranSIM(vector<double> valInitParmSIM,vector<double> valInitParam0){
    cout<<"tab1. Valeurs des Parametres pour SIMULATION, par defaut"<<endl;
    cout<<"NomParmSIM\t\t\tNouvelle Valeur\t\tValeur,par defaut\t\t\tSignification"<<endl;
    cout<<"1) tmax\t\t\t\t"<<valInitParmSIM[itmax]<<"\t\t\t"<<valInitParam0[itmax]<<"\t\t\t"<<"Temps pour fonctionner une simulation"<<endl;
    cout<<"2) nbVilles\t\t\t"<<valInitParmSIM[inbVille]<<"\t\t\t"<<valInitParam0[inbVille]<<"\t\t\t"<<"Nombre de villes"<<endl;
   // cout<<"3) rinitN\t\t\t"<<valInitParmSIM[irinitN]<<"\t\t\t"<<valInitParam0[irinitN]<<"\t\t\t"<<"Nombre de population initial"<<endl;
    //cout<<"4) rinitI\t\t\t"<<valInitParmSIM[irinitI]<<"\t\t\t"<<valInitParam0[irinitI]<<"\t\t\t"<<"Nombre de personnes infectées initial"<<endl;
   // cout<<"5) beta\t\t\t\t"<<valInitParmSIM[ibeta]<<"\t\t\t"<<valInitParam0[ibeta]<<"\t\t\t"<<"Taux de contact S->E initial"<<endl;
    cout<<"3) sigma\t\t\t\t"<<valInitParmSIM[isigma]<<"\t\t"<<valInitParam0[isigma]<<"\t\t"<<"Taux de transmission E->I"<<endl;
    cout<<"4) gamma\t\t\t\t"<<valInitParmSIM[igamma]<<"\t\t"<<valInitParam0[igamma]<<"\t\t"<<"Taux de transmission I->R"<<endl;
    cout<<"5) rmu\t\t\t\t"<<valInitParmSIM[irmu]<<"\t\t"<<valInitParam0[irmu]<<"\t\t"<<"Taux de mortalite et aussi le taux de connaissance"<<endl;
    cout<<"6) epsilon\t\t\t\t"<<valInitParmSIM[iepsilon]<<"\t\t\t"<<valInitParam0[iepsilon]<<"\t\t\t"<<"Taux de l'infection de l'exterieur"<<endl;
    cout<<"7) topology\t\t\t"<<valInitParmSIM[itopology]<<"\t\t\t"<<valInitParam0[itopology]<<"\t\t\t"<<"Type de topology entre des villes"<<endl;
    cout<<"8) rho\t\t\t"<<valInitParmSIM[irho]<<"\t\t\t"<<valInitParam0[irho]<<"\t\t\t"<<"Force de couplage entre deux villes"<<endl;
    cout<<"9) unitTemps\t\t\t"<<valInitParmSIM[iunitTemps]<<"\t\t\t"<<valInitParam0[iunitTemps]<<"\t\t\t"<<"Unit de temps pour sauvegarder les resultats"<<endl;
    cout<<"10) graine\t\t\t"<<valInitParmSIM[igraine]<<"\t\t\t"<<valInitParam0[igraine]<<"\t\t\t"<<"Chiffre de graine"<<endl;
    //iS0 = 10, iI0 = 11, iE0 = 12, iR0 = 13, iN0 = 14,
    //ibeta0 = 15, ibeta1 = 16, iphi0 = 17;
    cout<<"11) S0\t\t\t"<<valInitParmSIM[iS0]<<"\t\t\t"<<valInitParam0[iS0]<<"\t\t\t"<<"Nombre de susceptibles au moment t=0"<<endl;
    cout<<"12) E0\t\t\t"<<valInitParmSIM[iE0]<<"\t\t\t"<<valInitParam0[iE0]<<"\t\t\t"<<"Nombre d'exposés au moment t=0"<<endl;
    cout<<"13) I0\t\t\t"<<valInitParmSIM[iI0]<<"\t\t\t"<<valInitParam0[iI0]<<"\t\t\t"<<"Nombre d'infectés au moment t=0"<<endl;
    cout<<"14) R0\t\t\t"<<valInitParmSIM[iR0]<<"\t\t\t"<<valInitParam0[iR0]<<"\t\t\t"<<"Nombre de gueris au moment t=0"<<endl;
    cout<<"15) N0\t\t\t"<<valInitParmSIM[iN0]<<"\t\t\t"<<valInitParam0[iN0]<<"\t\t\t"<<"Nombre de population au moment t=0"<<endl;
    cout<<"16) beta0\t\t\t"<<valInitParmSIM[ibeta0]<<"\t\t\t"<<valInitParam0[ibeta0]<<"\t\t\t"<<"Valeur moyenne du taux de contact "<<endl;
    cout<<"17) beta1\t\t\t"<<valInitParmSIM[ibeta1]<<"\t\t\t"<<valInitParam0[ibeta1]<<"\t\t\t"<<"Amplitude du taux de contact"<<endl;
    cout<<"18) phiMIN\t\t\t"<<valInitParmSIM[iphiMIN]<<"\t\t\t"<<valInitParam0[iphiMIN]<<"\t\t\t"<<"Phase du forcage minimale"<<endl;
    cout<<"19) phiMAX\t\t\t"<<valInitParmSIM[iphiMAX]<<"\t\t\t"<<valInitParam0[iphiMAX]<<"\t\t\t"<<"Phase du forcage maximale"<<endl;



    }

double nouvelValSIM(int indexParam, double veilVal){
    double val;
    switch(indexParam)
    {
        case 0: cout<<"tmax = ";
                val = donnerChoix(); return val; break;
        case 1: cout<<"nbVille = ";  val = donnerChoix(); return val; break;
       // case 2: cout<<"rinitN = "; val = donnerChoix(); return val; break;
       // case 3: cout<<"rinitI = "; val = donnerChoix(); return val; break;
       // case 4: cout<<"beta = "; val = donnerChoix(); return val; break;
        case 2: cout<<"sigma = "; val = donnerChoix(); return val; break;
        case 3: cout<<"gamma = "; val = donnerChoix(); return val; break;
        case 4: cout<<"rmu = "; val = donnerChoix(); return val; break;
        case 5: cout<<"epsilon = "; val = donnerChoix(); return val; break;
        case 6: cout<<"topologie = "; val = donnerChoix(); return val; break;
        case 7: cout<<"rho = "; val = donnerChoix(); return val; break;
        case 8: cout<<"unitTemps = "; val = donnerChoix(); return val; break;
        case 9: cout<<"graine = "; val = donnerChoix(); return val; break;
        case 10: cout<<"S0 = "; val = donnerChoix(); return val; break;
        case 11: cout<<"E0 = "; val = donnerChoix(); return val; break;
        case 12: cout<<"I0 = "; val = donnerChoix(); return val; break;
        case 13: cout<<"R0 = "; val = donnerChoix(); return val; break;
        case 14: cout<<"N0 = "; val = donnerChoix(); return val; break;
        case 15: cout<<"beta0 = "; val = donnerChoix(); return val; break;
        case 16: cout<<"beta1 = "; val = donnerChoix(); return val; break;
        case 17: cout<<"phiMIN = "; val = donnerChoix(); return val; break;
        case 18: cout<<"phiMAX = "; val = donnerChoix(); return val; break;
        default : cout<<"Il n'y a aucune change!"<<endl; return veilVal; break;
      }
 }

vector<double> changerValInitParmSIM(vector<double> valInitParmSIM){
    vector<double> valInitParam0 = valInitParmSIM;
    coutEcranSIM(valInitParmSIM,valInitParam0);
    int choixChangeParam = opterChangerParam();
    if(choixChangeParam == 1){
        return valInitParmSIM;
    }
    else if(choixChangeParam == 2){
        coutEcranSIM(valInitParmSIM,valInitParam0);
        int flag = 1;
        while(flag == 1){
            cout<<"Quelle parametres que vous voulez changer, indexParam = ";
            int indexParam = donnerChoix();
            indexParam = indexParam - 1;
            valInitParmSIM[indexParam] = nouvelValSIM(indexParam,valInitParmSIM[indexParam]);
            coutEcranSIM(valInitParmSIM,valInitParam0);
            cout<<"Est-ce que vous voulez continuer à changer les valeurs des Parametres?(0/1), choix = ";
            flag = donnerChoix();
            if(flag != 1) {
                cout<<"Finir la change des valeurs de parametres!"<<endl;
                break;
            }
        }
        return valInitParmSIM;
    }
    else{
        cout<<"Programme marche avec les valeurs par defaut!"<<endl;
        return valInitParmSIM;
    }

//    return valInitParmSIM;
}

/******************************************************************************************************/
//Parametres de VACCINATION
//Initialiser les valeurs initiales pour les parametres
vector<double> initialervalParVAC(int argc, char** argv){
    vector<double> res;
    //itstartVAC = 0, itfinalVAC = 1, iPeriode = 2, itauxglobVAC = 3;
    double tstartVAC = 0.0;
    res.push_back(tstartVAC);// tmax = 2000(jours)
    double tfinalVAC = 2000.0;
    res.push_back(tfinalVAC); // nbVill = 1
    double Periode = 10.0;
    res.push_back(Periode);
    double tauxglobVAC = 80.0;
    res.push_back(tauxglobVAC);
    double nbMaxVaccins = 2000;
    res.push_back(nbMaxVaccins);

   //Parameters pour conserver les entrees
    char str[255] = "";
   for (int i=1; i<argc; i++)
       {
           // Recuperer le premier point de la fonction de transformation
             strcpy(str,argv[i]);
             if ((strcmp(str,"-tstartVAC") == 0))
             {
                 i++;
                 res[tstartVAC] =(double) atof(argv[i]);
             }
             else if((strcmp(str,"-tfinalVAC")) == 0)
             {
                 i++;
                 res[itfinalVAC] = atof(argv[i]);
             }
             else if((strcmp(str,"-Periode"))==0)
             {
                 i++;
                 res[iPeriode] = atof(argv[i]);
             }
             else if((strcmp(str,"-tauxglobVAC"))==0)
             {
                 i++;
                 res[itauxglobVAC] = atof(argv[i]);
             }
             else if((strcmp(str,"-nbMaxVaccins"))==0)
             {
                 i++;
                 res[inbMaxVaccins] = atof(argv[i]);
             }

        }

return res;

}

void coutEcranVAC(vector<double> valInitParmVAC,vector<double> valInitParam0){
    cout<<"tab2. Valeurs des Parametres pour VACCINATION, par defaut"<<endl;
    cout<<"NomParmVAC\t\t\tNouvelle Valeur\t\tValeur,par defaut\t\t\tSignification"<<endl;
    cout<<"1.tstartVAC\t\t\t\t"<<valInitParmVAC[itstartVAC]<<"\t\t\t"<<valInitParam0[itstartVAC]<<"\t\t\t"<<"Temps où on commencer la politique de vaccination"<<endl;
    cout<<"2.tfinalVAC\t\t\t\t"<<valInitParmVAC[itfinalVAC]<<"\t\t\t"<<valInitParam0[itfinalVAC]<<"\t\t\t"<<"Temps où on finir la politique de vaccination"<<endl;
    cout<<"3.Periode\t\t\t\t"<<valInitParmVAC[iPeriode]<<"\t\t\t"<<valInitParam0[iPeriode]<<"\t\t\t"<<"Periode entre deux fois de vaccination"<<endl;
    cout<<"4.tauxglobVAC\t\t\t\t"<<valInitParmVAC[itauxglobVAC]<<"\t\t\t"<<valInitParam0[itauxglobVAC]<<"\t\t\t"<<"Pourcentage de susceptibles vaccines"<<endl;
    cout<<"5.nbMaxVaccins\t\t\t\t"<<valInitParmVAC[inbMaxVaccins]<<"\t\t\t"<<valInitParam0[inbMaxVaccins]<<"\t\t\t"<<"Nombre de vaccins maximal"<<endl;

    }

double nouvelValVAC(int indexParam, double veilVal){
    double val;
    switch(indexParam)
    {
        case 0: cout<<"tstartVACC = "; val = donnerChoix(); return val; break;
        case 1: cout<<"tfinalVAC = ";  val = donnerChoix(); return val; break;
        case 2: cout<<"Periode = "; val = donnerChoix(); return val; break;
        case 3: cout<<"tauxglobVAC = "; val = donnerChoix(); return val; break;
        case 4: cout<<"nbMaxVaccins = "; val = donnerChoix(); return val; break;
        default : cout<<"Il n'y a aucune change!"<<endl; return veilVal; break;
      }
 }

vector<double> changerValInitParmVAC(vector<double> valInitParmVAC){
    vector<double> valInitParam0 = valInitParmVAC;
    coutEcranVAC(valInitParmVAC,valInitParam0);
    int choixChangeParam = opterChangerParam();
    if(choixChangeParam == 1){
        return valInitParmVAC;
    }
    else if(choixChangeParam == 2){
        coutEcranVAC(valInitParmVAC,valInitParam0);
        int flag = 1;
        while(flag == 1){
            cout<<"Quelle parametres que vous voulez changer, indexParam = ";
            int indexParam = donnerChoix();
            indexParam = indexParam - 1;
            valInitParmVAC[indexParam] = nouvelValVAC(indexParam,valInitParmVAC[indexParam]);
            coutEcranVAC(valInitParmVAC,valInitParam0);
            cout<<"Est-ce que vous voulez continuer à changer les valeurs des Parametres?(0/1), choix = ";
            flag = donnerChoix();
            if(flag != 1) {
                cout<<"Finir la change des valeurs de parametres!"<<endl;
                break;
            }
        }
        return valInitParmVAC;
    }
    else{
        cout<<"Programme marche avec les valeurs par defaut!"<<endl;
        return valInitParmVAC;
    }

//    return valInitParmSIM;
}

/******************************************************************************************************/
//Parametres de EXPORATIOn
//Initialiser les valeurs initiales pour les parametres
vector<double> initialervalParEVOL(int argc, char** argv){
    vector<double> res;
    //
    double nbEXPLOR = 5;
    res.push_back(nbEXPLOR);
    double nbSTRAT = 20;
    res.push_back(nbSTRAT);
    double nbPOLIT = 100;
    res.push_back(nbPOLIT);

   //Parameters pour conserver les entrees
    char str[255] = "";
   for (int i=1; i<argc; i++)
       {
           // Recuperer le premier point de la fonction de transformation
             strcpy(str,argv[i]);
             if ((strcmp(str,"-nbEXPLOR") == 0))
             {
                 i++;
                 res[inbEXPLOR] =(double) atof(argv[i]);
             }
             else if((strcmp(str,"-nbSTRAT")) == 0)
             {
                 i++;
                 res[inbSTRAT] = atof(argv[i]);
             }
             else if((strcmp(str,"-nbPOLIT"))==0)
             {
                 i++;
                 res[inbPOLIT] = atof(argv[i]);
             }

        }

return res;

}

void coutEcranEVOL(vector<double> valInitParmEVOL,vector<double> valInitParam0){
    cout<<"tab1. Valeurs des Parametres pour EVOLUTION, par defaut"<<endl;
    cout<<"NomParmEVOL\t\t\tNouvelle Valeur\t\tValeur,par defaut\t\t\tSignification"<<endl;
    cout<<"1.nbEXPLOR\t\t\t\t"<<valInitParmEVOL[inbEXPLOR]<<"\t\t\t"<<valInitParam0[inbEXPLOR]<<"\t\t\t"<<"Nombre de fois d'exporations"<<endl;
    cout<<"2.nbSTRAT\t\t\t\t"<<valInitParmEVOL[inbSTRAT]<<"\t\t\t"<<valInitParam0[inbSTRAT]<<"\t\t\t"<<"Nombre de strategies dans une fois d'exporations"<<endl;
    cout<<"3.nbPOLIT\t\t\t\t"<<valInitParmEVOL[inbPOLIT]<<"\t\t\t"<<valInitParam0[inbPOLIT]<<"\t\t\t"<<"Nombre de politiques essayees dans une strategie"<<endl;

    }

double nouvelValEVOL(int indexParam, double veilVal){
    double val;
    switch(indexParam)
    {
        case 0: cout<<"nbEXPLOR = ";val = donnerChoix(); return val; break;
        case 1: cout<<"nbSTRAT = "; val = donnerChoix(); return val; break;
        case 2: cout<<"nbPOLIT = "; val = donnerChoix(); return val; break;
        default : cout<<"Il n'y a aucune change!"<<endl; return veilVal; break;
      }
 }

vector<double> changerValInitParmEVOL(vector<double> valInitParmEVOL){
    vector<double> valInitParam0 = valInitParmEVOL;
    coutEcranEVOL(valInitParmEVOL,valInitParam0);
    int choixChangeParam = opterChangerParam();
    if(choixChangeParam == 1){
        return valInitParmEVOL;
    }
    else if(choixChangeParam == 2){
        coutEcranEVOL(valInitParmEVOL,valInitParam0);
        int flag = 1;
        while(flag == 1){
            cout<<"Quelle parametres que vous voulez changer, indexParam = ";
            int indexParam = donnerChoix();
            indexParam = indexParam - 1;
            valInitParmEVOL[indexParam] = nouvelValEVOL(indexParam,valInitParmEVOL[indexParam]);
            coutEcranEVOL(valInitParmEVOL,valInitParam0);
            cout<<"Est-ce que vous voulez continuer à changer les valeurs des Parametres?(0/1), choix = ";
            flag = donnerChoix();
            if(flag != 1) {
                cout<<"Finir la change des valeurs de parametres!"<<endl;
                break;
            }
        }
        return valInitParmEVOL;
    }
    else{
        cout<<"Programme marche avec les valeurs par defaut!"<<endl;
        return valInitParmEVOL;
    }

//    return valInitParmSIM;
}

/******************************************************************************************************/
//Le but de cette fonction est d'initialiser les nombres de prosonnes dans chaque groupe (S, E, I, R, N) pour chaque ville
//  nbVill: nombre de villes
//  neq: nombre de groupes
//  y: table pour sauvegarder les nombres de personnes dans chaque groupe
//  S0: nombre de suscptible t=0
//  E0: nombre d'exposes t=0
//  I0: nombre d'infecte t=0
//  R0: nombre de gueris t=0
//  N0: nombre de population t=0
void initialerY_00(int nbVill,int neq, unsigned long **&y, unsigned long S0, unsigned long E0, unsigned long I0,
                                                        unsigned long R0, unsigned long N0){
    y = new unsigned long*[nbVill];
    for (int nextville = 0; nextville <  nbVill; nextville++)
    {
        y[nextville] = new unsigned long[neq];
        y[nextville][iI] = I0;//rnitN*0.01;
        y[nextville][iN] = N0;
        y[nextville][iS] = S0;//y[nextville][iN]*0.2;//-y[nextville][iI];
        y[nextville][iE] = E0;
        y[nextville][iR] = R0;
        /*
        y[nextville][iI] = rnitI;//rnitN*0.01;
        y[nextville][iN] = rnitN;
        y[nextville][iS] = rnitN - rnitI;//y[nextville][iN]*0.2;//-y[nextville][iI];
        y[nextville][iE] = 0;
        y[nextville][iR] = 0;//y[nextville][iN]-y[nextville][iS]-y[nextville][iI];
        */

    }
    return;
}

void initialerY_01(int nbVill,int neq, unsigned long **&y, vector<nodeSEIR>*valeursSEIR0, unsigned long S0, unsigned long E0, unsigned long I0,
                                                        unsigned long R0, unsigned long N0){
    y = new unsigned long*[nbVill];
    for (int i = 0; i <  nbVill; i++)
    {
        y[i] = new unsigned long[neq];
        if(valeursSEIR0){
        y[i][iS] = (*valeursSEIR0)[i].S;
        y[i][iE] = (*valeursSEIR0)[i].E;
        y[i][iI] = (*valeursSEIR0)[i].I;
        y[i][iR] = (*valeursSEIR0)[i].R;
        y[i][iN] = (*valeursSEIR0)[i].N;
        }
        else{
            y[i][iI] = I0;//rnitN*0.01;
            y[i][iN] = N0;
            y[i][iS] = S0;//y[nextville][iN]*0.2;//-y[nextville][iI];
            y[i][iE] = E0;
            y[i][iR] = R0;

        }
    }
    return;
}
/******************************************************************************************************/
//Retourne un nombre pseudo-aléatoire selon une loi normale de paramètres mu et sigma
// @param mu moyenne (espérance mathématique) de la distribution
 // @param sigma écart-type de la distribution (doit être strictement positif)
double creerBetaAleatoire(double beta0, double sigma2) {

        // On récupère deux nombres pseudo-aléatoires indépendants selon une loi uniforme sur l'intervalle [0;1]
        double randNumUni = ((double) rand())/((double) RAND_MAX);
        double randNumBi = ((double) rand())/((double) RAND_MAX);

        // On récupère un nombre pseudo-aléatoire selon une loi normale centrée réduite
        // (Paramètres : moyenne = 0, écart-type = 1)
        // Utilisation de l'algorithme de Box-Muller
        double cos2Pi = cos(DEUX_PI*randNumBi);
        double randNumNorm = 0.0;
        if(cos2Pi>= 0.0)
        {
            randNumNorm = sqrt(-2.0*log(randNumUni))*cos2Pi;
        }
        else{
            randNumNorm = sqrt(-2.0*log(randNumUni))*(-cos2Pi);
        }
        return (beta0 + sigma2 * randNumNorm);
}

/******************************************************************************************************/
// Fonction à calculer le taux de contact beta_i qui subit un forcage saisonier, pour la population i.
// beta_i(t) = beta0 *(1+ beta1*cos(2*pi*t _ phi_i);
double* calculerBeta_01(int nbVill, double beta0, double beta1, double t, double *phi){
    double *res = new double[nbVill];
    for(int i=0; i<nbVill; i++){
        res[i] = 0.0;
        res[i] = beta0*(1 + beta1 *cos(2*Pi*t + phi[i]));
    }
    return res;
}

/******************************************************************************************************/
// Fonction à calculer les phases pour chaque la population
double *calculerPhases_01(int nbVill, double phiMAX, double phiMIN){
    double *res = new double[nbVill];
    double unit_phi = 0.0;
    if(nbVill==1) unit_phi = phiMIN;
    else
     unit_phi = (double)(phiMAX-phiMIN)/(nbVill-1);
    double phi_tp=0.0;
    for(int i=0; i<nbVill; i++){
        res[i] = phi_tp;
        phi_tp = phi_tp + unit_phi;
        //cout<<res[i]<<endl;
    }

    return res;
}

/******************************************************************************************************/
// Fonction à calculer les valeurs de couplage pour entre deux population i et j
double **calculerTauxCouplage_01(int nbVill, double rho0){
    double **arr_rho;
    arr_rho = new double*[nbVill];
    for (int i = 0; i<nbVill; i++)
    {
        arr_rho[i] = new double[nbVill];
        for(int j = 0; j< nbVill ; j++) {
            if(j==i) arr_rho[i][j] = 1.0;
            else
            arr_rho[i][j] = rho0;
           // cout<<arr_rho[i][j]<<endl;
        }
    }
    return arr_rho;
}

/******************************************************************************************************/
vector<nodeSEIR> readTable(string filename){
    vector<nodeSEIR> data;
    FILE *pFile;
    char mystring [255];
    nodeSEIR nodeSEIRTemps;
    double t;
    int S, E, I,R,N;
    stringstream ss ;
    pFile = fopen (filename.c_str(), "r");
    if (pFile == NULL) perror ("Error opening file in the readTable function!");
    else {
        while (!feof(pFile)) {
        fgets (mystring , 255 , pFile);
        ss.clear ();
        ss.str ("");
        ss<<mystring;
        ss>>t >>S >>E >>I >>R>>N;
       // cout<<"t = "<<t<<" S = "<<S<<" E = "<<E<< " I="<<I<<endl;
        nodeSEIRTemps.t=t;
        nodeSEIRTemps.S = S;
        nodeSEIRTemps.E = E;
        nodeSEIRTemps.I = I;
        nodeSEIRTemps.R = R;
        nodeSEIRTemps.N = N;
        nodeSEIRTemps.event = 0;
        data.push_back(nodeSEIRTemps);
    }
        fclose (pFile);
    }
    return data;
}

/******************************************************************************************************/
// Extraire les series SEIR[x, x+T] dans la partie sationnaire de I
vector<nodeSEIR> extraireSEIR(double x, double T, vector<nodeSEIR> data)
{
    // Executer la fonctionner 1.4
    // Extraire SEIR[x, x+T] dans la parite stationnaire de I
    double xExtrait = x;
    vector<nodeSEIR> dataExtrait;
    double xT = xExtrait + T;
    int sizeData = data.size();

    // Creer un fichier pour sauvegarder les series SEIR[x, x+T]
    FILE * fSEIR;
    char nomfSEIR[255];
    sprintf(nomfSEIR,"Output/SEIRx%.1fT%.1f.csv",xExtrait, T);
    fSEIR = fopen (nomfSEIR,"w+");
    char lSEIR[255];

    // Prendre les series SEIR [x, x+T]
    nodeSEIR nodeTpSEIR;
    for(int i=0; i<sizeData; i++)
    {
        if((data[i].t >= xExtrait)&(data[i].t <= xT)){
            nodeTpSEIR.t = data[i].t - xExtrait;
            nodeTpSEIR.S = data[i].S;
            nodeTpSEIR.E = data[i].E;
            nodeTpSEIR.I = data[i].I;
            nodeTpSEIR.R = data[i].R;
            nodeTpSEIR.N = data[i].N;
            nodeTpSEIR.event = data[i].event;
            dataExtrait.push_back(nodeTpSEIR);
            sprintf(lSEIR, "%.1f         %lu    %lu     %lu     %lu    %lu\n",data[i].t , data[i].S, data[i].E,  data[i].I, data[i].R,data[i].N);
            fputs (lSEIR,fSEIR);
        }
    }

   fclose (fSEIR);
   return dataExtrait;

}

/******************************************************************************************************/
// Extraire les valeurs initiales de SEIR pour les villes dans les series SEIR selon la valeur de synchronie ySYNCH
vector<nodeSEIR> extraireSEIR0(double ySYNCH,double tstartSYNCH,double T,int nbVill,
                               vector<nodeSEIR> dataExtrait){
    vector<nodeSEIR> vectSEIR0;
    // Sauvegarder les valeurs initiales des S, E, I, R, N des populations dans le fichier "initSEIR....txt"
    FILE * fileInitialSEIR;
    char nomfInit[255];
    sprintf(nomfInit,"Output/initSEIRy%.3fT%.1ftstart%.2f_SYNCH.csv",ySYNCH,T,tstartSYNCH);
    fileInitialSEIR = fopen (nomfInit,"w+");
    char lineInitialSEIR[255];
    fputs ("TempsParJour   S     E      I     R      N\n",fileInitialSEIR);
    cout<<"ysynchronie= "<<ySYNCH<<  "   periode= "<<T<<endl;
    double    tstep = (double)(ySYNCH*T/nbVill);
    cout<<"tstep = "<<tstep<<endl;
    int szData = dataExtrait.size();
    tstartSYNCH = dataExtrait[0].t;
    for(int i=0; i< nbVill; i++){
        double point = tstartSYNCH + i*tstep;
        int pos=0;
        for(int j=0; j<szData; j++)
        {
            if(dataExtrait[j].t >= point) {
                pos = j;
                break;
            }
        }
        vectSEIR0.push_back(dataExtrait[pos]);
        sprintf(lineInitialSEIR, "%.5f         %lu    %lu     %lu     %lu    %lu\n",point ,dataExtrait[pos].S, dataExtrait[pos].E,dataExtrait[pos].I,  dataExtrait[pos].R,dataExtrait[pos].N);
        cout<<"ville "<<i<< "  "<< lineInitialSEIR<<endl;
        fputs (lineInitialSEIR,fileInitialSEIR);
    }
    fclose(fileInitialSEIR);
    return vectSEIR0;
}

/******************************************************************************************************/
// Vaccination
// but est de vacciner p% de la population S->R
// Avec tous les T jours,
// La même chose pour toutes les villes
int vaccinerToutesVilles_V1(double tauxglobVAC,int nbVilles, unsigned long **y,double_vector *politique)
{
    unsigned long nbPersonneVaccine=0;
        int i;
cout<<"In vaccinerToutesVilles_V1"<<endl;
    if (politique == NULL){
        cout<<"politique == NULL"<<endl;
                for(int i=0;i<nbVilles;i++){
                        nbPersonneVaccine=(y[i][iS]*tauxglobVAC)/100;
                        y[i][iS] = y[i][iS]-nbPersonneVaccine;
                        y[i][iR] = y[i][iR]+nbPersonneVaccine;
                }
    }
        else
    {
        cout<<"politique != NULL"<<endl;
                for(i=0;i<nbVilles;i++)
                        if (y[i][iS]*(*politique)[0]+y[i][iE]*(*politique)[1]+y[i][iI]*(*politique)[2]+y[i][iR]*(*politique)[3]+100.0*(*politique)[4] > 0) {
//			if (y[i][iS]*(*politique)[0]+100.0*(*politique)[1] > 0) {
                                        nbPersonneVaccine=(y[i][iS]*tauxglobVAC)/100;
                                        y[i][iS] = y[i][iS]-nbPersonneVaccine;
                                        y[i][iR] = y[i][iR]+nbPersonneVaccine;
                        }
    }

        return nbPersonneVaccine;
}

/******************************************************************************************************/
// C'est la fonction principale simulerSEIR_00
vector<string> simulerSEIR_00(vector<double> valParSIM, char SAUVER[5],
                           vector<nodeSEIR>*valeursSEIR0,
                           vector<double>*valParVAC,vector<double>*politiqueEVOL)
{
    // pour SIMULATION
    // Initialer le nombre de cas = 5 (S, E, I, R, N); le nombre d'evenements = 9
    int neq = 6, nevent = 9;
    // le nombre de villes utilisé
    int nbVill = valParSIM[inbVille];
    // Type de Topology :
    // topology = 0, modèle en ^iles
    // topology = 1, modèle en cercle
    // topology = 2 , modèle en réseaux
    int topology = valParSIM[itopology];
    //Taux de E->I
    double sigma =valParSIM[isigma];
    //Taux de I->R
    double gamma = valParSIM[igamma];
    //Taux de transmission
    double beta = valParSIM[ibeta0];
    //Taux de natalité
    double rmu = valParSIM[irmu];
    //Taux de transmission entre deux villes
    double rho = valParSIM[irho];
    double epsilon = valParSIM[iepsilon];
    //Parametres de probabilités, p1: probabilité de temps, p2: probabilité de ville, p3: probabilité de réaction qui se produit
    double p1 = 0.0, p2 = 0.0, p3 = 0.0;
    double tstep = 0.0;
    // Initialisation les nombres de personnes dans le groupe I et le groupe N
    unsigned long rinitN = valParSIM[iN0];
    unsigned long rinitI = valParSIM[iI0];
    // Initialisation
    double tmax = valParSIM[itmax];
    double fsum = 0.0, fsumVilles = 0.0;
    double t = 0.0;
    int m = 0, nextville = 0;
    double fprop = 0.0, fvill = 0.0, fvillTemps = 0.0;
    double unitTemps = valParSIM[iunitTemps];
    double graine = valParSIM[igraine];
    int sizeTable = 0, sizeTblDer = 0;
    //int MAXSIZEEVENTTABLE = 5000;
    unsigned long S0 = valParSIM[iS0];
    unsigned long E0 = valParSIM[iE0];
    unsigned long I0 = valParSIM[iI0];
    unsigned long R0 = valParSIM[iR0];
    unsigned long N0 = valParSIM[iN0];

    // pour  VACCINATION
    int periode = 0.0;
    double tempsPolitiqueVaccination = 0.0;
    int kMax = 0;
    //Le temps où on commence à vacciner
    double tstart = 0.0;
    //Le temps où on finit à vaccinaer
    double tfinal = 0.0;
    int tpNombreFrequence = 0;
    double pointVacciner = 0.0;
    int nombreSusceptibleVaccinees = 0;
    //double tauxVaccToutesVilles = (*valParVAC)[itauxglobVAC];
    //int nombreMaxVaccins = (*valParVAC)[inbMaxVaccins];
   // int nombreTotalVaccine = 0;
   // int nombreTotalInfectes = 0;
    if(valParVAC){
        periode = (*valParVAC)[iPeriode];
        tempsPolitiqueVaccination = 0.0;
        kMax = 0;
        //Le temps où on commence à vacciner
        tstart = (*valParVAC)[itstartVAC];
        //Le temps où on finit à vaccinaer
        tfinal = (*valParVAC)[itfinalVAC];
        tpNombreFrequence = 0;
        pointVacciner = tstart + tpNombreFrequence*periode;
        nombreSusceptibleVaccinees = 0;
        //tauxVaccToutesVilles = 10;
        tauxVaccToutesVilles = (*valParVAC)[itauxglobVAC];
        nombreMaxVaccins = (*valParVAC)[inbMaxVaccins];
        nombreTotalVaccine = 0;
        nombreTotalInfectes = 0;
        tempsPolitiqueVaccination = tfinal-tstart;
        kMax = (int)(tempsPolitiqueVaccination /periode) +1;
    }


    // La programme marche :
    //Initialisation
    double *lamda = new double[nbVill];

    //Sauvegarder les valeurs dans chaque cas et dans chaque evenement
    unsigned long **y;
    double **f;
    f = new double*[nbVill];
    for (int nextville = 0; nextville <  nbVill; nextville++)
    {
        f[nextville] = new double[nevent];
        for(int i = 0; i< nevent ; i++) f[nextville][i] = 0.0;
    }


    //Creer les tables qui sauvegardent des resultats quand le programme marche
    unsigned long table[nbVill][5000][7];
    unsigned long indexTable[nbVill];
    //
    unsigned long tableSauverDernierLigne[nbVill][5000][7];
    unsigned long indexTblDer[nbVill];

    for(int i=0; i<nbVill;i++){
        indexTable[i]=0;
        indexTblDer[i]=0;

    }

    //Creer une variable de temps
    double pointTemps = unitTemps;
    //Creer les FILES
    string listeNomFichier[nbVill];
    FILE *totalFile;
    char nameFileTotal[255];
    char line0[255] ="t(jours)\t S \t E \t I \t R \t N \t event \n";
    char lineTotal0[255] ="t(jours)\t totalS \t totalE \t totalI \t totalR \n";

    if(valeursSEIR0 == NULL)//Simulateur SEIR sans synchronie
       {
           initialerY_00(nbVill, neq, y, S0,E0, I0, R0, N0);
          if(valParVAC == NULL){
              //Initialiser les fichiers pour sauvegarder ler resultal par ville
              //ville i: t    S    E   I   R    N
              //        0.001 10   20  20  20   70
              //        ...    ..  ..  ..  ..   ..
              int nbNomVilles=0;
              if((strcmp(SAUVER,"Y"))==0){
                  char nameFile[255];
                  FILE * pFile;
                  for(int i = 0; i < nbVill;i++){
                      int ville = i+1;
                      sprintf(nameFile,"Output/BatchVil%d_NbVil%dEps%.3fTop%d_SIM.csv",ville,nbVill,rho,topology);
                      pFile = fopen (nameFile,"w+");
                      listeNomFichier[nbNomVilles]= string(nameFile);
                      fputs(line0,pFile);
                      fclose(pFile);
                      nbNomVilles++;
                  }
              }
              //Creer un nom de fichier total qui sauvegarde le total de S, E, I, R au moment t
              sprintf(nameFileTotal,"Output/BatchtolSEIRNbVil%dEps%.3fTop%dGam%.2fSig%.2f_SIM.csv",nbVill,rho,topology,gamma,sigma);
          }
          else
          {
              //Initialiser les fichiers pour sauvegarder ler resultal par ville
              //ville i: t    S    E   I   R    N
              //        0.001 10   20  20  20   70
              //        ...    ..  ..  ..  ..   ..
              int nbNomVilles=0;
              if((strcmp(SAUVER,"Y"))==0){
                  char nameFile[255];
                  FILE * pFile;
                  for(int i = 0; i < nbVill;i++){
                      int ville = i+1;
                      sprintf(nameFile,"Output/BatchVil%d_NbVil%dEps%.3fTop%d_VAC.csv",ville,nbVill,rho,topology);
                      pFile = fopen (nameFile,"w+");
                      listeNomFichier[nbNomVilles]= string(nameFile);
                      fputs(line0,pFile);
                      fclose(pFile);
                      nbNomVilles++;
                  }
              }
              //Creer un nom de fichier total qui sauvegarde le total de S, E, I, R au moment t
              sprintf(nameFileTotal,"Output/BatchtolSEIRNbVil%dEps%.3fTop%dGam%.2fSig%.2f_VAC.csv",nbVill,rho,topology,gamma,sigma);

          }


       }
        else //Simulateur avec la synchronie
        {
            //vector<nodeSEIR> SEIR0 = *valeursSEIR0;
            y = new unsigned long*[nbVill];
            for (int i = 0; i <  nbVill; i++)
            {
                y[i] = new unsigned long[neq];
                y[i][iS] = (*valeursSEIR0)[i].S;
                y[i][iE] = (*valeursSEIR0)[i].E;
                y[i][iI] = (*valeursSEIR0)[i].I;
                y[i][iR] = (*valeursSEIR0)[i].R;
                y[i][iN] = (*valeursSEIR0)[i].N;
            }

            //Initialiser les fichiers pour sauvegarder ler resultal par ville
            //ville i: t    S    E   I   R    N
            //        0.001 10   20  20  20   70
            //        ...    ..  ..  ..  ..   ..
            int nbNomVilles=0;
            if((strcmp(SAUVER,"Y"))==0){
                char nameFile[255];
                FILE * pFile;
                for(int i = 0; i < nbVill;i++){
                    int ville = i+1;
                    sprintf(nameFile,"Output/BatchVil%d_nbVil%dEps%.3fTop%d_SYNCH.csv",ville,nbVill,rho,topology);
                    pFile = fopen (nameFile,"w+");
                    listeNomFichier[nbNomVilles]= string(nameFile);
                    fputs(line0,pFile);
                    fclose(pFile);
                    nbNomVilles++;
                }
            }

            //Creer un fichier total qui sauvegarde le total de S, E, I, R au moment t
            sprintf(nameFileTotal,"Output/BatchtolSEIRNbVil%dEps%.3fTop%dGam%.2fSig%.2f_SYNCH.csv",nbVill,rho,topology,gamma,sigma);
        }
    totalFile = fopen(nameFileTotal,"w+");
    fputs(lineTotal0,totalFile);

    //Sauvegarder les valeurs des parametres dans un fichier BatchParmVil....txt pour chaque ville
    for(int i = 0; i < nbVill; i++){
        streambuf *psbufTOTAL, *backupTOTAL;
        ofstream filestrTOTAL;
        char nomBatchParmTOTAL[255];
        sprintf(nomBatchParmTOTAL,"Output/BatchParmVil%d_nbVil%dEps%.3fTop%d.txt",i+1,nbVill,rho,topology);
        filestrTOTAL.open (nomBatchParmTOTAL);
        backupTOTAL = cout.rdbuf();
        psbufTOTAL = filestrTOTAL.rdbuf();
        cout.rdbuf(psbufTOTAL);
        cout<< " Parametres pour SIMULATION"<<endl;
        cout << " tmax \t\t" << tmax<<endl;
        cout << " nbVill \t"<< nbVill<<endl;
        cout << " rinitS \t"<<y[i][iS]<<endl;
        cout << " rinitE \t"<<y[i][iE]<<endl;
        cout << " rinitI \t"<<y[i][iI]<<endl;
        cout << " rinitR \t"<<y[i][iR]<<endl;
        cout << " rinitN \t"<<y[i][iN]<<endl;
        cout << " beta \t\t"<<beta<<endl;
        cout << " gamma \t\t"<<gamma<<endl;
        cout << " sigma \t\t"<<sigma<<endl;
        cout << " rmu \t\t"<<rmu<<endl;
        cout << " epsilon \t\t"<<epsilon<<endl;
        cout << " topologie \t "<<topology<<endl;
        cout << " rho \t"<<rho<<endl;
        cout << " unitTemps \t"<<unitTemps<<endl;
        cout << " graine \t"<<graine<<endl;
        if(valParVAC){
            cout << " Parametres pour VACCINATION"<<endl;
            cout << " tstartVAC \t"<<(*valParVAC)[itstartVAC]<<endl;
            cout << " tfinalVAC \t"<<(*valParVAC)[itfinalVAC]<<endl;
            cout << " Periode \t"<<(*valParVAC)[iPeriode]<<endl;
            cout << " tauxglobVAC \t"<<(*valParVAC)[itauxglobVAC]<<endl;
            cout << " nbMaxVaccins \t"<<(*valParVAC)[inbMaxVaccins]<<endl;
        }
        //else
          //  cout << " Programme ne fait pas la vaccination"<<endl;
        cout.rdbuf(backupTOTAL);
        filestrTOTAL.close();
        //
    }


    //Calculer le temps où la programme marche selon le temps de CPU
    clock_t start, end;
    start = clock();

    //Initialiser les nombres de personnes de chaque groupe
    //srand(graine);
    //initialerY( nbVill, neq, y, rinitN,  rinitI);

    //Initialiser les valeurs des tables
    // à la fois sauvegarder les valeurs initiale des parametres et variables pour chaque ville
     for (int i = 0; i <  nbVill; i++)
     {
            table[i][0][it] = 0;
            table[i][0][iS] = y[i][iS];
            table[i][0][iE] = y[i][iE];
            table[i][0][iI] = y[i][iI];
            table[i][0][iR] = y[i][iR];;
            table[i][0][iN] = y[i][iN];
            table[i][0][ievent] = 0;
            indexTable[i]++;

            tableSauverDernierLigne[i][0][it] =0;
            tableSauverDernierLigne[i][0][iS] = y[i][iS];
            tableSauverDernierLigne[i][0][iE] = y[i][iE];
            tableSauverDernierLigne[i][0][iI] = y[i][iI];
            tableSauverDernierLigne[i][0][iR] = y[i][iR];
            tableSauverDernierLigne[i][0][iN] = y[i][iN];
            tableSauverDernierLigne[i][0][ievent] = 0;
            indexTblDer[i]++;

     }

    //Initialer les valeurs de Lamda pour chaque ville
    calculerLamda_00(nbVill, topology,beta,rho,y,lamda);

    //Initialiser les valeurs pour la fonction de propensité
    calculerF(nbVill,epsilon,rmu,sigma,gamma,lamda,y,f);

    //Caculer le temps où la réaction se produit
    while(t<tmax){

        //Generate uniform random numbers
        p1 = rand()/(double)RAND_MAX;
        p2 = rand()/(double)RAND_MAX;
        p3 = rand()/(double)RAND_MAX;

        //Caculer le total fsum de toutes les villes
        fsumVilles = calculerFsumToutesVilles(nbVill,nevent,f);

        //Determine time interval and update time
        if(fsumVilles > 0.0) {
            tstep = -log(p1/fsumVilles)/fsumVilles;
            //tstep = -log(p1)/fsumVilles;
            t = t + tstep;
        }
        else {
            cout<<"no event"<<endl;
            cout<<"parce que: fsumVilles = "<<fsumVilles<<endl;
            break;
        }

        //Faire vacciner ici si c'est le temps de vaccination
        //Vacciner pour toutes les villes selon chaque la taux de vaccination par ville
        if(valParVAC){
            if((t>pointVacciner)&&(tpNombreFrequence<= kMax) && nombreTotalVaccine < nombreMaxVaccins){
                //cout<<"t="<<t<< " pointVacciner="<<pointVacciner<<" kMax="<<kMax<<endl;
                nombreSusceptibleVaccinees=vaccinerToutesVilles_V1(tauxVaccToutesVilles,nbVill,y,politiqueEVOL);
                nombreTotalVaccine += nombreSusceptibleVaccinees;
                pointVacciner += periode;
                tpNombreFrequence++;
            }
        }
        //Selectionne aleatoirement une ville ou se produira l'evenement
        fvillTemps =0;
        for(int i=0; i<nbVill; i++){
            fvill = calculerFsum(nevent,i,f);
            fvillTemps = fvillTemps+fvill;
            if(p2 < fvillTemps/fsumVilles) {
                nextville = i;
                break;
            }
        }


        // Selectionne aleatoirement la nature du prochain evenement dans nextville
        fprop = 0.0;
        m = -1;
        fsum = calculerFsum(nevent,nextville,f);
        for(int j=0; j< nevent; j++){
            fprop = fprop+f[nextville][j];
            if(p3 < (fprop/fsum)) {
                m = j;
                break;
            }
        }

        // Faire l'événement m
        if(m !=-1){
                //Faire l'événement m
                fairEvenementM(m,nextville,y);

                //Recalculer les valeurs de Lamda
                //beta = creerBetaAleatoire(beta0,sigma2);
                calculerLamda_00(nbVill,topology,beta,rho,y,lamda);

                //Recalculer les valleurs de la foction de propensité
                calculerF(nbVill,epsilon,rmu,sigma,gamma,lamda,y,f);

                //Recalculer la valeur de fsum
                fsum= calculerFsum(nevent,nextville,f);

                //Caculer le total fsum de toutes les villes
                fsumVilles = calculerFsumToutesVilles(nbVill,nevent,f);

                //Sauvegarder la valeur de t, le nombre de personnes dans le groupe I et l'événement qui se produit
                if(t> pointTemps){
                    sizeTable = indexTable[nextville];

                    if(sizeTable < 5000){
                        table[nextville][sizeTable][0] = pointTemps;
                        table[nextville][sizeTable][1] = y[nextville][iS];
                        table[nextville][sizeTable][2] = y[nextville][iE];
                        table[nextville][sizeTable][3] = y[nextville][iI];
                        table[nextville][sizeTable][4] = y[nextville][iR];
                        table[nextville][sizeTable][5] = y[nextville][iN];
                        table[nextville][sizeTable][6] = m;
                        indexTable[nextville]++;
                    }
                    else{
                        if((strcmp(SAUVER,"Y"))==0) {
                            writeInFile(listeNomFichier[nextville],nextville,table,indexTable);
                        }
                        sizeTblDer = indexTblDer[nextville];
                        for(int i=0; i<7; i++){
                            tableSauverDernierLigne[nextville][sizeTblDer][i]= table[nextville][sizeTable][i];
                        }
                        indexTable[nextville]=0;
                    }

                    pointTemps +=unitTemps;

                    //Sauvegarder dans le fichier total
                    unsigned long totalS=0, totalI=0, totalE=0, totalR=0;
                    char line[255];
                    for(int k=0; k<nbVill; k++){
                        unsigned long dernierObject = indexTable[k];
                        if(dernierObject>0){
                            dernierObject = dernierObject-1;
                            totalS += table[k][dernierObject][iS];
                            totalE += table[k][dernierObject][iE];
                            totalI += table[k][dernierObject][iI];
                            totalR += table[k][dernierObject][iR];
                    }
                        else{
                            dernierObject = indexTblDer[k]-1;
                            totalS += table[k][dernierObject][iS];
                            totalE += table[k][dernierObject][iE];
                            totalI += table[k][dernierObject][iI];
                            totalR += table[k][dernierObject][iR];
                        }
                    }
                    sprintf(line, " %.0f \t\t %ld \t\t %ld \t\t %ld \t\t %ld \n",t, totalS, totalE, totalI, totalR);
                    fputs (line,totalFile);
                }


            }

        else
        {
             // PB !!
               cout<<"Giang a pas reussi car il y a un m=-1!"<<endl;
        }

        //Le programme fini!
    }

    end = clock();

    // Le temps CPU utilisé:
    double totalDure = (double)(end - start)/CLOCKS_PER_SEC;
    //cout<<"Le temps CPU utilisé est de :  "<<totalDure<<" s"<<endl;

    //Afficher les resultats
    for(int nextville=0; nextville<nbVill; nextville++){
        sizeTable = indexTable[nextville];
        sizeTblDer = indexTblDer[nextville];
        for(int i=0; i<7; i++){
            tableSauverDernierLigne[nextville][sizeTblDer][i]= table[nextville][sizeTable][i];
        }
    }
   // AfficheDerniereValeur(nbVill,rinitN,rho,tmax,tableSauverDernierLigne,indexTblDer);


    //Afficher tous les resultats de chaque ville
    if((strcmp(SAUVER,"Y"))==0){
        for(int nextville=0; nextville<nbVill; nextville++){
            writeInFile(listeNomFichier[nextville],nextville,table,indexTable);
        }
    }

    //Afficher tous les arguments de l'appel
    //afficherArgument(nbVill,rinitN,rho,tmax,argc, argv);

    //Libérer la mémoire
    delete []lamda;
    for (int nextville = 0; nextville <  nbVill; nextville++)
    {
        delete []y[nextville];
        delete []f[nextville];

    }
    delete []y;
    delete[]f;
   // delete[]table;
    fclose(totalFile);
    //cout<<"Le programme est avec succes!"<<endl;
    vector<string> arNomFichierRes;
        for(int i=0; i<nbVill;i++)
        {
            //cout<<"name i=  "<<i <<"  "<<listeNomFichier[i]<<endl;
            arNomFichierRes.push_back((string)listeNomFichier[i].c_str());
        }
        arNomFichierRes.push_back((string)nameFileTotal);
        return(arNomFichierRes);
    //The end the programme principal
    }
/******************************************************************************************************/
string *creerListeNomFichier(vector<double> valParSIM){
    // le nombre de villes utilisé
    int nbVill = valParSIM[inbVille];
    // Type de Topology :
    // topology = 0, modèle en ^iles
    // topology = 1, modèle en cercle
    // topology = 2 , modèle en réseaux
    int topology = valParSIM[itopology];
    //Taux de E->I
    double sigma =valParSIM[isigma];
    double tmax = valParSIM[itmax];
    //Taux de I->R
    double gamma = valParSIM[igamma];
    //Taux de transmission
    double beta0 = valParSIM[ibeta0];
    double beta1 = valParSIM[ibeta1];
    //Taux de natalité
    double rmu = valParSIM[irmu];
    //Taux de transmission entre deux villes
    double epsilon = valParSIM[iepsilon];
    double phi0 = valParSIM[iphiMIN];

    double unitTemps = valParSIM[iunitTemps];
    double graine = valParSIM[igraine];

    unsigned long S0 = valParSIM[iS0];
    unsigned long E0 = valParSIM[iE0];
    unsigned long I0 = valParSIM[iI0];
    unsigned long R0 = valParSIM[iR0];
    unsigned long N0 = valParSIM[iN0];

    //Taux de la force de couplage entre deux villes i et j
    //double rho = valParSIM[irho];
    double rho0 = valParSIM[irho];

    string listeNomFichier[nbVill];
    char line0[255] ="t(jours)\t S \t E \t I \t R \t N \t event \n";
    //Initialiser les fichiers pour sauvegarder ler resultal par ville
    //ville i: t    S    E   I   R    N
    //        0.001 10   20  20  20   70
    //        ...    ..  ..  ..  ..   ..
    int nbTp=0;
        char nameFile[255];
        FILE * pFile;
        for(int i = 0; i < nbVill;i++){
            int ville = i+1;
            sprintf(nameFile,"Output/SimVil%d_NbVil%dTmax%.3fSigma%.3fGamma%.3fRmu%.3fEpsilon%.3fTop%dRho%.3fUnittemps%.3fGraine%.3fS%dlE%dI%dlR%dN%dlbeta0_%.3fbeta1_%.3fphi0_%.3f.csv",ville,nbVill,tmax,sigma,gamma,rmu,epsilon,topology,rho0,unitTemps,graine,S0,E0,I0,R0,N0,beta0,beta1,phi0);
            pFile = fopen (nameFile,"w+");
            listeNomFichier[nbTp]= string(nameFile);
            fputs(line0,pFile);
            fclose(pFile);
            nbTp++;
        }

   return listeNomFichier;
}

/******************************************************************************************************/
// C'est la fonction principale simulerSEIR_01
vector<string> simulerSEIR_01(vector<double> valParSIM, char SAUVER[5],
                           vector<nodeSEIR>*valeursSEIR0, double *phi, double **arr_rho){
    // pour SIMULATION
    // Initialer le nombre de cas = 5 (S, E, I, R, N); le nombre d'evenements = 9
    int neq = 6, nevent = 9;
    // le nombre de villes utilisé
    int nbVill = valParSIM[inbVille];
    // Type de Topology :
    // topology = 0, modèle en ^iles
    // topology = 1, modèle en cercle
    // topology = 2 , modèle en réseaux
    int topology = valParSIM[itopology];
    //Taux de E->I
    double sigma =valParSIM[isigma];
    //Taux de I->R
    double gamma = valParSIM[igamma];
    //Taux de transmission
    double beta0 = valParSIM[ibeta0];
    double beta1 = valParSIM[ibeta1];
    //Taux de natalité
    double rmu = valParSIM[irmu];

    //Taux de transmission entre deux villes
    double epsilon = valParSIM[iepsilon];
    double phiMIN = valParSIM[iphiMIN];
    double phiMAX = valParSIM[iphiMAX];

    unsigned long S0 = valParSIM[iS0];
    unsigned long E0 = valParSIM[iE0];
    unsigned long I0 = valParSIM[iI0];
    unsigned long R0 = valParSIM[iR0];
    unsigned long N0 = valParSIM[iN0];

    //Taux de la force de couplage entre deux villes i et j
    //double rho = valParSIM[irho];
    double rho0 = valParSIM[irho];


    //Parametres de probabilités, p1: probabilité de temps, p2: probabilité de ville, p3: probabilité de réaction qui se produit
    double p1 = 0.0, p2 = 0.0, p3 = 0.0;
    double tstep = 0.0;

    // Initialisation
    double tmax = valParSIM[itmax];
    double fsum_i = 0.0, fsum_global = 0.0;
    double t = 0.0;
    int evenement = 0, nextville = 0;
    double fprop = 0.0, fvill = 0.0, fvillTemps = 0.0;
    double unitTemps = valParSIM[iunitTemps];
    double graine = valParSIM[igraine];
    int sizeTable = 0, sizeTblDer = 0;
    //int MAXSIZEEVENTTABLE = 5000;

    cout<<"tmax = "<<tmax<<", nbVille = "<<nbVill<<", sigma = "<<sigma<<", gamma = "<<gamma<<", rmu = "<<rmu<<", epsilon = "<<epsilon<<endl;
    cout<<"topology = "<<topology<<", rho0 = "<<rho0 <<", unitTemps = "<<unitTemps<<", graine = "<<graine<<endl;
    cout<<"S0 = "<<S0<<"I0 = "<<I0<<"E0 = "<<E0<<", R0 = "<<R0<<", N0 = "<<N0<<endl;
    cout<<"beta0 = "<<beta0<<", beta1 = "<<beta1<<", phiMIN = "<<phiMIN<<", phiMAX = "<<phiMAX<< endl;


    // La programme marche :
    //Initialisation
    double *lamda = new double[nbVill];
    double *beta = new double[nbVill];
   // double *phi  = new double[nbVill];
    for(int i = 0; i< nbVill ; i++) {
        lamda[i] = 0.0;
        beta[i] = 0.0;
      //  phi[i] = 0.0;
    }


    //Sauvegarder les valeurs dans chaque cas et dans chaque evenement
    unsigned long **y;
    initialerY_01(nbVill,neq,y,valeursSEIR0,S0,E0,I0,R0,N0);
  //  for(int i=0; i<nbVill; i++){
    //    cout<<"ville i="<<i<< " S= "<<y[i][iS]<<"  E="<<y[i][iE]<<" I="<<y[i][iI]<<" R="<<y[i][iR]<<"  N="<<y[i][iN]<<endl;
    //}


    //Sauvegarder les valeur de fonctions de propabilite
    double **f;
    f = new double*[nbVill];
    for (int nextville = 0; nextville <  nbVill; nextville++)
    {
        f[nextville] = new double[nevent];
        for(int i = 0; i< nevent ; i++) f[nextville][i] = 0.0;
    }

    //Creer les tables qui sauvegardent des resultats quand le programme marche
    unsigned long table[nbVill][5000][7];
    unsigned long indexTable[nbVill];
    //
    unsigned long tableSauverDernierLigne[nbVill][5000][7];
    unsigned long indexTblDer[nbVill];

    for(int i=0; i<nbVill;i++){
        indexTable[i]=0;
        indexTblDer[i]=0;

    }


    //Creer une variable de temps
    double pointTemps = unitTemps;

    //Creer les FILES
    string listeNomFichier[nbVill];
    if((strcmp(SAUVER,"Y"))==0){
        char line0[255] ="t(jours)\t S \t E \t I \t R \t N \t event \n";
        //Initialiser les fichiers pour sauvegarder ler resultal par ville
        //ville i: t    S    E   I   R    N
        //        0.001 10   20  20  20   70
        //        ...    ..  ..  ..  ..   ..
        int nbTp=0;
            char nameFile[255];
            FILE * pFile;
            for(int i = 0; i < nbVill;i++){
                int ville = i+1;
                sprintf(nameFile,"Output/SimVil%d_NbVil%dTmax%.3fSigma%.3fGamma%.3fRmu%.3fEpsilon%.3fTop%dRho%.3fUnittemps%.3fGraine%.3fS%ldE%ldI%ldR%ldN%dbeta0_%.3fbeta1_%.3fphiMIN%.3fphiMAX%.3f.csv",ville,nbVill,tmax,sigma,gamma,rmu,epsilon,topology,rho0,unitTemps,graine,S0,E0,I0,R0,N0,beta0,beta1,phiMIN,phiMAX);
                pFile = fopen (nameFile,"w+");
                listeNomFichier[nbTp]= string(nameFile);
                fputs(line0,pFile);
                fclose(pFile);
                nbTp++;
            }
    }

    //for(int i = 0; i< nbVill ; i++) {
      //  cout<<"nom= "<<i<<listeNomFichier[i] <<endl;
    //}
    //Creer un nom de fichier total qui sauvegarde le total de S, E, I, R au moment t
    FILE *totalFile;
    char nameFileTotal[255];
    char lineTotal0[255] ="t(jours)\t totalS \t totalE \t totalI \t totalR \n";
    sprintf(nameFileTotal,"Output/BatchtolSEIRNbVil%dTmax%.3fSigma%.3fGamma%.3fRmu%.3fEpsilon%.3fTop%dRho%.3fUnittemps%.3fGraine%.3fS%ldE%ldI%ldR%ldN%ldbeta0_%.3fbeta1_%.3fphiMIN%.3fphiMAX%.3f.csv",nbVill,tmax,sigma,gamma,rmu,epsilon,topology,rho0,unitTemps,graine,S0,E0,I0,R0,N0,beta0,beta1,phiMIN,phiMAX);
    totalFile = fopen(nameFileTotal,"w+");
    fputs(lineTotal0,totalFile);

    //Calculer le temps où la programme marche selon le temps de CPU
    clock_t start, end;
    start = clock();

    //Initialiser les nombres de personnes de chaque groupe

   // srand(graine);

    //Initialiser les valeurs des tables
    // à la fois sauvegarder les valeurs initiale des parametres et variables pour chaque ville
     for (int i = 0; i <  nbVill; i++)
     {
            table[i][0][it] = 0;
            table[i][0][iS] = y[i][iS];
            table[i][0][iE] = y[i][iE];
            table[i][0][iI] = y[i][iI];
            table[i][0][iR] = y[i][iR];;
            table[i][0][iN] = y[i][iN];
            table[i][0][ievent] = 0;
            indexTable[i]++;

            tableSauverDernierLigne[i][0][it] =0;
            tableSauverDernierLigne[i][0][iS] = y[i][iS];
            tableSauverDernierLigne[i][0][iE] = y[i][iE];
            tableSauverDernierLigne[i][0][iI] = y[i][iI];
            tableSauverDernierLigne[i][0][iR] = y[i][iR];
            tableSauverDernierLigne[i][0][iN] = y[i][iN];
            tableSauverDernierLigne[i][0][ievent] = 0;
            indexTblDer[i]++;

     }

    //Caculer le temps où la réaction se produit
    while(t<tmax){

        //Initialer les valeurs de beta pour chaque ville
         beta = calculerBeta_01(nbVill,beta0,beta1,t,phi);

        //Calculer les valeurs de Lamda
         calculerLamda_01(nbVill,epsilon,beta,arr_rho,y,lamda);
         /*
         for(int i = 0; i< nbVill ; i++) {
             cout<<"i= "<<i<<lamda[i] <<endl;
              cout<<"i= "<<i<<beta[i] <<endl;
           cout<<"i= "<<i<<phi[i] <<endl;
         }
         */

        //Calculer les valleurs de la foction de propensité
        calculerF(nbVill,epsilon,rmu,sigma,gamma,lamda,y,f);

        //Generate uniform random numbers
        p1 = rand()/(double)RAND_MAX;
        p2 = rand()/(double)RAND_MAX;
        p3 = rand()/(double)RAND_MAX;

        //Caculer le total fsum de toutes les villes
        fsum_global = calculerFsumToutesVilles(nbVill,nevent,f);

        //Determine time interval and update time
        if(fsum_global > 0.0) {
            tstep = -log(p1)/fsum_global;
            t = t + tstep;
        }
        else {
            cout<<"no event"<<endl;
            cout<<"parce que: fsumVilles = "<<fsum_global<<endl;
            break;
        }

        //Selectionne aleatoirement une ville ou se produira l'evenement
        fvillTemps =0;
        for(int i=0; i<nbVill; i++){
            fvill = calculerFsum(nevent,i,f);
            fvillTemps = fvillTemps+fvill;
            if(p2 <= fvillTemps/fsum_global) {
                nextville = i;
                break;
            }
        }

        // Selectionne aleatoirement la nature du prochain evenement dans nextville
        fprop = 0.0;
        evenement = -1;
        fsum_i = calculerFsum(nevent,nextville,f);
        for(int j=0; j< nevent; j++){
            fprop = fprop+f[nextville][j];
            if(p3 <= (fprop/fsum_i)) {
                evenement = j;
                break;
            }
        }

        // Faire l'événement m
        if(evenement !=-1){
                //Faire l'événement m
                fairEvenementM(evenement,nextville,y);

                //Sauvegarder la valeur de t, le nombre de personnes dans le groupe I et l'événement qui se produit
                if(t> pointTemps){
                    sizeTable = indexTable[nextville];
                    if(sizeTable < 5000){
                        table[nextville][sizeTable][0] = pointTemps;
                        table[nextville][sizeTable][1] = y[nextville][iS];
                        table[nextville][sizeTable][2] = y[nextville][iE];
                        table[nextville][sizeTable][3] = y[nextville][iI];
                        table[nextville][sizeTable][4] = y[nextville][iR];
                        table[nextville][sizeTable][5] = y[nextville][iN];
                        table[nextville][sizeTable][6] = evenement;
                        indexTable[nextville]++;
                    }
                    else{
                        if((strcmp(SAUVER,"Y"))==0) {
                            writeInFile(listeNomFichier[nextville],nextville,table,indexTable);
                        }
                        sizeTblDer = indexTblDer[nextville];
                        for(int i=0; i<7; i++){
                            tableSauverDernierLigne[nextville][sizeTblDer][i]= table[nextville][sizeTable][i];
                        }
                        indexTable[nextville]=0;
                    }

                    pointTemps +=unitTemps;

                    //Sauvegarder dans le fichier total
                    unsigned long totalS=0, totalI=0, totalE=0, totalR=0;
                    char line[255];
                    for(int k=0; k<nbVill; k++){
                        unsigned long dernierObject = indexTable[k];
                        if(dernierObject>0){
                            dernierObject = dernierObject-1;
                            totalS += table[k][dernierObject][iS];
                            totalE += table[k][dernierObject][iE];
                            totalI += table[k][dernierObject][iI];
                            totalR += table[k][dernierObject][iR];
                    }
                        else{
                            dernierObject = indexTblDer[k]-1;
                            totalS += table[k][dernierObject][iS];
                            totalE += table[k][dernierObject][iE];
                            totalI += table[k][dernierObject][iI];
                            totalR += table[k][dernierObject][iR];
                        }
                    }
                    sprintf(line, " %.0f \t\t %ld \t\t %ld \t\t %ld \t\t %ld \n",t, totalS, totalE, totalI, totalR);
                    fputs (line,totalFile);
                }


            }

        else
        {
             // PB !!
               cout<<"Giang a pas reussi car il y a un m=-1!"<<endl;
        }

        //Le programme fini!
    }

    end = clock();

    // Le temps CPU utilisé:
    double totalDure = (double)(end - start)/CLOCKS_PER_SEC;
    //cout<<"Le temps CPU utilisé est de :  "<<totalDure<<" s"<<endl;

    //Afficher les resultats
    for(int nextville=0; nextville<nbVill; nextville++){
        sizeTable = indexTable[nextville];
        sizeTblDer = indexTblDer[nextville];
        for(int i=0; i<7; i++){
            tableSauverDernierLigne[nextville][sizeTblDer][i]= table[nextville][sizeTable][i];
        }
    }
   // AfficheDerniereValeur(nbVill,rinitN,rho,tmax,tableSauverDernierLigne,indexTblDer);


    //Afficher tous les resultats de chaque ville
    if((strcmp(SAUVER,"Y"))==0){
        for(int nextville=0; nextville<nbVill; nextville++){
            writeInFile(listeNomFichier[nextville],nextville,table,indexTable);
        }
    }

    //Afficher tous les arguments de l'appel
    //afficherArgument(nbVill,rinitN,rho,tmax,argc, argv);

    //Libérer la mémoire
    delete []lamda;
    for (int nextville = 0; nextville <  nbVill; nextville++)
    {
        delete []y[nextville];
        delete []f[nextville];

    }
    delete []y;
    delete[]f;
   // delete[]table;
    fclose(totalFile);
    //cout<<"Le programme est avec succes!"<<endl;
    vector<string> arNomFichierRes;
        for(int i=0; i<nbVill;i++)
        {
            //cout<<"name i=  "<<i <<"  "<<listeNomFichier[i]<<endl;
            arNomFichierRes.push_back((string)listeNomFichier[i].c_str());
        }
        arNomFichierRes.push_back((string)nameFileTotal);
        return(arNomFichierRes);
    //The end the programme principal
}
/******************************************************************************************************/
//Le but de cette fonction est de calculer les valeurs de Lamda (c'est le taux de transmission S->E) pour chaque ville
//  nbVill: nombre de villes
//  topology: type de topology
//  beta: valeur de trasmission
//  rho: taux de transmission entre deux villes
//  y: table pour sauvegarder les nombres de personnes dans chaque groupe
//  lamda: table pour sauvegarder les valeurs de lamda
void calculerLamda_00(int nbVill, int topology, double beta, double rho,
                   unsigned long **y, double *lamda){
    double sumIj;
    double tp,tpi;
    double tp1,tp2;

    for(int v=0; v<nbVill; v++){
        if(topology == 0){// en forme en ^iles
            sumIj = 0;
        }

        if(topology==1) {//en forme de cycle
            sumIj=0.0;
            int villefinal = nbVill-1;
            if(villefinal==0){//il y a seulement une ville
                sumIj = y[v][iI];
            }
            else if(villefinal==1)//il y a deux villes
            {
                if(v==0) sumIj=y[v+1][iI];
                else sumIj =y[v-1][iI];
            }
            else //il y a plus de deux villes
            {
            if(v==0){
                tp1 = (double)y[v+1][iI];
                tp2 = (double)y[villefinal][iI];
                sumIj= sumIj+tp1+tp2;
            }
            else if(v==villefinal){
                tp1 = (double)y[v-1][iI];
                tp2 = (double)y[0][iI];
                sumIj=sumIj+tp1+tp2;
            }
            else{
                tp1 = (double)y[v-1][iI];
                tp2 = (double)y[v+1][iI];
                sumIj=sumIj+tp1+tp2;
            }
          }
        }

        if(topology ==2){//en forme d'un graphe complet
            sumIj=0.0;//
            for(int l=0; l<nbVill;l++){
                if(l !=v){
                    tp = (double)y[l][iI];
                    sumIj= sumIj+tp;
                }
            }
        }

        tpi =(double)y[v][iI];
        lamda[v] = (double)beta*(tpi+ rho*sumIj)/y[v][iN];
    }
}

/******************************************************************************************************/
//Le but de cette fonction est de calculer les valeurs de Lamda (c'est le taux de transmission S->E) pour chaque ville
//  nbVill: nombre de villes
//  beta: valeur de trasmission
//  rho: taux de transmission entre deux villes i et j
//  y: table pour sauvegarder les nombres de personnes dans chaque groupe
//  lamda: table pour sauvegarder les valeurs de lamda
//  epsilon: taux d'infection de l'exterieur
void calculerLamda_01(int nbVill, double epsilon, double *beta, double **rho,
                   unsigned long **y, double *lamda){
    double total_ij = 0.0;
    for(int i=0; i<nbVill; i++){
        total_ij = 0.0;
        for(int j=0; j<nbVill;j++){
            total_ij = (rho[i][j]*(epsilon*beta[i]*y[j][iN] + (1-epsilon)*beta[j]*y[i][iN])*y[j][iI])/(y[i][iN]*y[j][iN]);
        }
        lamda[i] = (double)(total_ij/nbVill);
    }

}

/******************************************************************************************************/
//Le but de cette fonction est de calculer les valeurs de la fonction de propensité
//  nbVill: nombre de villes
//  rmu: taux de naissance et taux de mortalité
//  sigma: taux de transmission E->I
//  gamma: taux de transmission I->R
//  lamda: table à sauvegarder les valeur de lamda
//  y:  table à  sauvegarder les nombres de personnes dans chaque groupe
//  f: table pour sauvegarder les valeurs de la fonction de propensité
void calculerF(int nbVill,double epsilon,double rmu, double sigma, double gamma, double *lamda,
               unsigned long **y,double **f){
    for(int nextville=0; nextville< nbVill;nextville++){

        //infection from outside
        //f[nextville][0]= epsilon;//*y[nextville][iS];


        //Birth
        f[nextville][1] = rmu*y[nextville][iN];

        //Susceptible death
        f[nextville][2] = rmu*y[nextville][iS];

        //Infection
        f[nextville][3] = lamda[nextville]*y[nextville][iS] + epsilon*y[nextville][iS];

        //Exposed death
        f[nextville][4] = rmu*y[nextville][iE];

        //Movement into next class
        f[nextville][5] = sigma*y[nextville][iE];

        //Infected death
        f[nextville][6] = rmu*y[nextville][iI];

        //Movement into next class
        f[nextville][7] = gamma*y[nextville][iI];

        //Recovered death
        f[nextville][8] = rmu*y[nextville][iR];
    }
}

/******************************************************************************************************/
//Le but de cette fonction est de calculer le total des fonctions de propensités pour une ville
//  nevent: nombre d'événements
//  nextville: ville i
//  f: table àsauvegarder les valeurs des fonctions de propensités
//  nevent: nombre d'événements
double calculerFsum(int nevent, int ville,double **f){
    double fsum=0.0;
    for(int i=0; i< nevent;i++){
        fsum=fsum + f[ville][i];
    }
    return fsum;
}

/******************************************************************************************************/
//Le but de cette fonction est de calculer le somme fsum de toutes les villes
double calculerFsumToutesVilles(int nbVill, int nevent, double **f){
    double fsumVilles =0.0;
    for(int i = 0; i<nbVill; i++){
        for(int j = 0; j<nevent; j++){
            fsumVilles = fsumVilles + f[i][j];
        }
    }
    return fsumVilles;
}

/******************************************************************************************************/
//Le but de cette fonction est de faire l'événement m
//  m: événement m
//  nextville: ville nextville
//  y: table à  sauvegarder les nombres de personnes dans chaque groupe
void fairEvenementM(int m, int nextville, unsigned long **y){

    //m=0 : infection de dehors
    /*

    if((m ==0)&&(y[nextville][iS]>0)){
         y[nextville][iI]=y[nextville][iI]+1;
         y[nextville][iS]=y[nextville][iS]-1;
         nombreTotalInfectes++;

    }
    */


    // événement : une personne est née
    if(m ==1){
        y[nextville][iS]=y[nextville][iS]+1;
        y[nextville][iN]=y[nextville][iN]+1;
    }

    // événement : une personne est morte
    if(m ==2){
        y[nextville][iS] = y[nextville][iS]-1;
        y[nextville][iN]=y[nextville][iN]-1;
    }

    // événement : une personne est exposée S->E
    if(m == 3){
        y[nextville][iS]=y[nextville][iS]-1;
        y[nextville][iE]=y[nextville][iE]+1;
    }

    // événement : une personne exposée est morte
    if(m==4){
        y[nextville][iE]=y[nextville][iE]-1;
        y[nextville][iN]=y[nextville][iN]-1;
    }

    // événement : une personne exposée est infecté E->I
    if(m == 5){
        y[nextville][iE]=y[nextville][iE]-1;
        y[nextville][iI]=y[nextville][iI]+1;
        nombreTotalInfectes++;
    }

    // événement : une personne infectée est morte
    if(m == 6){
        y[nextville][iI]=y[nextville][iI]-1;
        y[nextville][iN]=y[nextville][iN]-1;
    }

    // événement : une personne infecté est guérie
    if(m == 7){
        y[nextville][iI]=y[nextville][iI]-1;
        y[nextville][iR]=y[nextville][iR]+1;
    }

    // événement : une personne guérie est morte
    if(m==8){
        y[nextville][iR]=y[nextville][iR]-1;
        y[nextville][iN]=y[nextville][iN]-1;
    }
}

/******************************************************************************************************/
// Le but de cette fonction est de sauvegarder les resultats dans les fichier de text, un fichier par ville
// nameFile: le nom de fichier qui sauvegarde les resultats
// iVill: ville i
// table: table contient les resultat
 void writeInFile(string nameFile,int ivill,
                  unsigned long table[100][5000][7],unsigned long indexTable[100]){

    int sizeTable;
    char line[200];
    unsigned long S, E, I, R, N;
    int event;
    int t;
    FILE * pFile;
    sizeTable = indexTable[ivill];
    //int ville = ivill+1;
    pFile = fopen (nameFile.c_str(),"a+");
    for(int j=0; j<sizeTable;j++){

        t = table[ivill][j][0];
        S = table[ivill][j][1];
        E = table[ivill][j][2];
        I = table[ivill][j][3];
        R = table[ivill][j][4];
        N = table[ivill][j][5];
        event = table[ivill][j][6];
        sprintf(line, "%d \t\t %ld \t %ld \t %ld \t %ld \t %ld \t %d\n",t ,S, E, I, R, N, event);
        fputs (line,pFile);
    }
        fclose (pFile);
    return;
}

/******************************************************************************************************/
// Le but de cette fonction est de sauvegarder les dernieres lignes des resultats par ville
// Le nom de fichier a la forme suivante :
//    SimuNbVille<nombre des villes>Nini<nombre de population initial>
//        Epsion<valeur de rho>Tmax<valeur de tmax>_SUMMARY.csv
// nbVill: nombre de villes
// rinitI: le nombre de personnes infectees initiale
// rho: le taux de contact entre deux villes
// tmax: le temps maximal ou la simulation marche
// table: table contient les resultat
 void afficherDerniereValeur(int nbVill,unsigned long rinitN, double rho, double tmax,
                             unsigned long table[100][5000][7],unsigned long indexTabDer[100]){
    unsigned long S, E, I, R, N;//event;
    int t;
    unsigned long sizeTable;
    char line[200];
    FILE * pFile;
    int ville;
    unsigned long SIZE;
    char nameFile[100];
    sprintf(nameFile,"%s%d%s%ld%s%.3f%s%.1f%s","Output/SimuNbVille",nbVill,"Ninit",rinitN,"rho",rho,"Tmax",tmax,"_SUMMARY.csv");
    pFile = fopen (nameFile,"w");
    sprintf(line,"%s", "Villes;  S;    E;    I;     R;   N;   t\n");
    fputs (line,pFile);
    for(int ivill=0; ivill<nbVill; ivill++){
        SIZE = indexTabDer[ivill];
        sizeTable = SIZE-1;
        ville = ivill+1;
        S = table[ivill][sizeTable][iS];
        E = table[ivill][sizeTable][iE];
        I = table[ivill][sizeTable][iI];
        R = table[ivill][sizeTable][iR];
        N = table[ivill][sizeTable][iN];
        t = table[ivill][sizeTable][0];
        sprintf(line, "%d;      %ld;    %ld;     %ld;     %ld;     %ld;   %d  \n",ville,S, E, I, R, N,t);
        fputs (line,pFile);
    }
    fclose (pFile);
    return;
}

/******************************************************************************************************/
// Le but de cette fonction cree un fichier qui contient les valeurs de tous les paramètres en argument de l'appel
// Le nom de fichier a la forme suivante :
//          SimuNbVille<nombre des villes>Nini<nombre de population initial>
//                Epsion<valeur de rho>Tmax<valeur de tmax>_PARAMS.csv
// Le contenu du fichier a la forme suivante:
//              Arguments;  Nom;        Valeur
//              1;          Nbville;    1
//              2;          Nini;       5000
//              .......
void afficherArgument(int nbVill,unsigned long rinitN, double rho, double tmax,
                      int argc, char **argv){
    FILE *pFile;
    char nameFile[100];
    char line[200];
    sprintf(nameFile,"%s%d%s%ld%s%.3f%s%.1f%s","SimuNbVille",nbVill,"Ninit",rinitN,"rho",rho,"Tmax",tmax,"_PARAMS.csv");
    pFile = fopen (nameFile,"w");
    sprintf(line,"%s", "Argument;  Nom;    Valeur\n");
    fputs (line,pFile);

    int nbArgument=0;
    char str[50]="";
    for (int i=1; i<argc; i++)
    {
        // Recuperer le premier point de la fonction de transformation
          strcpy(str,argv[i]);
          nbArgument++;
          if ((strcmp(str,"-sigma") == 0))
          {
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"sigma", argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-gamma")) == 0)
          {
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"gamma",argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-beta"))==0)
          {
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"beta", argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-rmu"))==0)
          {
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"rmu", argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-initI"))==0)
          {
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"rInitI", argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-initN"))==0)
          {
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"rInitN", argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str, "-nbVill"))==0)
          {
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"nbVilles", argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-topology"))==0)
          {
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"topology", argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-rho"))==0){
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"rho", argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-tmax"))==0){
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"tmax", argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-unitTemps"))==0){
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"unitTemps",argv[i]);
              fputs (line,pFile);

          }
          else if((strcmp(str,"-SAUVER"))==0){
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"SAUVER",argv[i]);
              fputs (line,pFile);
          }
          else if((strcmp(str,"-graine"))==0){
              i++;
              sprintf(line,"%d;  %s;  %s\n",nbArgument,"graine",argv[i]);
              fputs (line,pFile);
          }
        }
    fclose(pFile);
    return;
}

/******************************************************************************************************/

// La but de cette fonction est de faire la fonction evo_start
void verifierEvolVac(vector<double> valParSIM, vector<double> valParVAC,
                     int nbEXPOL = 2,int nbSTRAT = 5,int nbPOL = 100){

    //Creer un fichier txt pour sauvegarder les resulats
    FILE * fEvolVac;
    char nomfEvolVac[255];
    sprintf(nomfEvolVac,"Output/ResEvolVacnbEXPOR%dnbSTRAT%d.csv",nbEXPOL, nbSTRAT);
    fEvolVac = fopen (nomfEvolVac,"w+");
    char lfEvolVac[255];
    sprintf(lfEvolVac, "Avantd'Optimisation \t\tApresd'Optimisation \n");
    fputs (lfEvolVac,fEvolVac);

    char SAUVER[5];
    SAUVER[0]='N';

    // mon objectif: reussir a diminuer le nombre d'infectes en utilisant 2000 vaccins au plus.
    // evaluate policy
    double_vector politiqueAPRES(5),fitness_sequenceAVANT, fitness_sequenceAPRES;
    for (int i=0; i < nbEXPOL; i++) {
            fitness_sequenceAVANT.clear();
            fitness_sequenceAPRES.clear();
            evo_strat(NULL,valParSIM,SAUVER,NULL,&valParVAC,nbSTRAT,1.0,true,&fitness_sequenceAVANT);
            politiqueAPRES = evo_strat(&politiqueAPRES,valParSIM,SAUVER,NULL,&valParVAC,nbSTRAT,1.0,true,&fitness_sequenceAPRES);

            for (int j=0; j < fitness_sequenceAPRES.size(); j++)
            {
                sprintf(lfEvolVac,"%.2f \t\t\t %.2f\n",fitness_sequenceAVANT[j],fitness_sequenceAPRES[j]);
                fputs (lfEvolVac,fEvolVac);
            }
    }
    fclose(fEvolVac);
}

/******************************************************************************************************/
// La but de cette fonction est de faire 100 fois d'essayer la fonction simulerSEIR
double eval_policy(vector<double> valParSIM,
                   char SAUVER[5],
                   vector<nodeSEIR>*valeursSEIR0,
                   vector<double>*valParVAC,
                   double_vector* policy)
{
        int trial;
        //int ntrials=100;
        int ntrials=100;
        int ev=0;
        for (trial=0; trial < ntrials; trial++) {
                //-- Make mymain be silent
                std::stringstream   redirectStream;
                std::streambuf*     old_cout  = std::cout.rdbuf( redirectStream.rdbuf() );
                simulerSEIR_00(valParSIM,SAUVER,valeursSEIR0,valParVAC,policy);
                //--- start talking again
                std::cout.rdbuf(old_cout);

                ev += nombreTotalInfectes;
        }

        //return - ev / ntrials;
        return - ev / ntrials;
}
/******************************************************************************************************/
